subroutine initpsi(x0,gama,x,psi,NSTATES,imr)
! Gera uma Função de Onda inicial na forma de uma gaussiana
implicit none
double precision, parameter :: nx = 10001
integer i, is, NSTATES, imr
double precision norm, gama, E, Ry, a0, am, ak, x0, dx, am1, am2, ec, pi
double complex :: img, dt
double precision, dimension(nx) :: x
double complex, dimension(nx,NSTATES) :: psi

	common / constants / dx, dt, img, Ry, a0, am, am1, am2, ec, pi

E = 0.D0/Ry !wavepacket energy in meV (divided by Ry to become dimensionless)

ak = dsqrt(am1*E) !wave vector (it is already dimensionless)

! *** Valor da Função de Onda Inicial
do is=1,NSTATES
	do i=1,nx
		psi(i,is) = cdexp(img*ak*x(i))*dexp(-(x(i)-x0)**2/(2.D0*(gama)**2))*x(i)**(is-1)
	enddo
enddo

! Normalização do Pacote Inicial
do is=1,NSTATES
	norm = 0.D0
	do i=1,nx
		norm = norm + cdabs(psi(i,is))**2*dx
	enddo
	do i=1,nx
		psi(i,is) = psi(i,is)/dsqrt(norm)
	enddo
enddo

if(imr.eq.1) then
	call GRAMSCHMIDT(psi,NSTATES)
end if

return
end subroutine initpsi
